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INTRODUCTION 

The  American  Cancer  Society  estimated  that  there  would  be  approximately  230,110  new  cases 
diagnosed  and  approximately  29,900  prostate  cancer  related  deaths  in  2004  (American  Cancer 
Society  2004)  Prostate  cancer  screening  today  generally  uses  the  Prostate  Specific  Antigen 
(PSA)  blood  testing,  free  PSA  testing  and  Digital  Rectal  Examination  (DRE).  When  using  a 
‘cutoff’  of  PSA  >  4.0  ng/mL  and  an  abnormal  DRE,  sensitivity,  specificity  and  Positive  Predictive 
Value  (PPV)  are  38%,  88%  and  56%  respectively  (Crawford,  Leewansangtong  et  al.  1999). 
When  either  an  elevated  PSA  or  an  abnormal  DRE  are  used,  (in  isolation  -  not  in  combination), 
sensitivity,  specificity  and  PPV  are  even  lower  (Crawford,  Leewansangtong  et  al.  1999).  When 
the  PSA  is  used  there  exists  a  significant  gray  area  (4-10  ng/mL)  in  which  cancers  may  be 
missed  and  yet  the  number  of  negative  biopsies  is  large.  Even  though  cancer  detection 
sensitivity,  specificity  and  PPV  are  improved  by  combining  PSA  and  DRE  (Toubert,  Schlageter 
et  al.  1990;  Crawford,  Leewansangtong  et  al.  1999)  the  usefulness  of  DRE  remains 
fundamentally  limited  due  to  its  subjective  nature.  Additionally,  DRE  is  practically  limited  to  the 
detection  of  shallow  (subcapsular)  palpable  abnormalities.  Even  systematic  multi-core  biopsy 
fails  to  detect  clinically  detectable  cancers  in  up  to  34%  of  men  (Frauscher,  Klauser  et  al.  2001). 
However,  there  is  evidence  that  as  additional  biopsies  cores  are  added,  sensitivity  improves 
(Taylor,  Gancarczyk  et  al.  2002).  This  observation  has  resulted  in  an  increase  the  number  of 
cores  taken  during  routine  examination.  Nevertheless,  biopsy-based  detection  sensitivity 
remains  less  than  ideal.  Thus,  there  is  plenty  of  compelling  clinical  interest  in  finding  improved 
methods  for  the  early  diagnosis  of  prostate  cancer  with  improved  sensitivity  and  specificity.  One 
recent  example  of  progress  in  the  field  of  prostate  cancer  detection  involves  an  effort  to 
automate  the  DRE  examination.  Savazyan  recently  described  a  system  for  'mechanical  imaging' 
of  the  prostate  (Sarvazyan  1998).  This  system  comprises  a  rectal  probe  that  is  instrumented 
with  an  array  of  pressure  sensing  strain  gages  and  a  3D  magnetic  positioner  device.  In  an  in 
vitro  trial  (Weiss,  Hartanto  et  al.  2001),  the  new  system  correctly  detected  and  located  100%  of 
the  nodules  under  examination.  This  compares  with  detection  rates  of  83%  and  67%  for  an 
experienced  urologist  and  a  student  respectively.  Thus,  a  significant  improvement  over  the 
conventional  DRE  examination  has  been  demonstrated  for  the  in  vitro  case.  Another  recent 
development  is  the  observation  that  the  sensitivity  of  an  ultrasound  examination  can  be 
improved  by  the  use  of  a  microbubble  based  contrast  agents  (Frauscher,  Klauser  et  al.  2001). 
Frauscher's  approach  (Frauscher,  Klauser  et  al.  2001)  involved  the  use  of  contrast  agent 
enhanced  Color  Doppler  that  improved  the  detection  of  hypervascular  regions  associated 
cancer.  Prostate  cancer  was  detected  by  contrast  agent  assisted  ultrasound  in  23  of  24  patients 
known  to  have  prostate  cancer.  (The  method  used  for  determining  definitively  which  patients 
had  cancer  is  not  entirely  clear  in  the  article.)  In  comparison,  conventional  ultrasound  detected 
cancer  in  17  patients.  The  contrast  agent  assisted  approach  detected  cancer  in  8  patients  with  a 
negative  systematic  biopsy-based  diagnosis.  However,  the  cost  of  the  contrast  agent  used  in 
this  study  was  $65  per  patient.  This  cost  makes  up  approximately  half  of  the  cost  of  a 
conventional  ultrasound  examination  and  therefore  represents  a  considerable  impediment  to  its 
widespread  acceptance.  However,  more  recent  publications  (Clements  2002;  Halpern,  McCue 
et  al.  2002)  (including  one  from  Frauscher's  group)  cast  doubt  on  the  true  extent  of  the 
improvement  in  diagnostic  accuracy  obtained  by  using  contrast  agents.  Specifically,  Halpern 
was  unable  to  detect  cancers  in  the  inner  gland  and  achieved  a  cancer  detection  sensitivity  of 
only  42%  (Halpern,  McCue  et  al.  2002). 
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BODY 

The  work  conducted  as  part  of  this  Army  funded  program  can  be  considered  as  divided  into  the  following 
key  “Aims”: 

1 .  Research,  design,  development  and  prototype  testing  of  a  new  transrectal  ultrasound  transducer, 
syringe  pump  and  ultrasound  instrumentation  to  facilitate  a  Synthetic  Digital  Rectal  Examination 
(SDRE). 

2.  Research,  development  and  prototype  testing  of  techniques  to  enable  quantitative  (dimensionally 
accurate)  3D  reconstructions  of  the  prostate 

3.  Research,  development  and  test  of  techniques  to  improve  ultrasound  image  quality  and  to 
facilitate  automated  (or  semi-automated)  border  detection  of  lesions 

4.  Small  scale  clinical  test  at  the  University  of  Virginia 

Progress  has  been  made  in  each  of  these  tasks  in  the  third  year  of  the  grant. 

Progress  with  respect  to  the  areas  are  related  directly  to  the  committed  Statement  of  Work  that  was 
funded: 

Aim  1.  Design,  specify,  and  have  built,  a  high  resolution  transducer  optimized  for 
imaging  elastic  inhomogeneities,  unsurpassed  B-Mode  image  resolution  and  possessing 
integrated  3D  capability. 

A  high  frequency  (8-14  MHz)  transducer  array  was  designed  and  specified  in  Year  1.  In  2005, 
the  transducer  was  manufactured  and  delivered  by  Vermon  SA,  Tours,  France.  This  specialized 
ultrasound  transducer  has  two  tracking  arrays  each  with  32  elements,  a  central  imaging  array 
with  192  elements  and  the  elements  are  spaced  on  a  0.2  mm  pitch.  This  transducer  is  providing 
superb  imaging  resolution  in  a  transducer  housing  specially  designed  for  transrectal  ultrasound. 
The  transducer  has  provided  the  high  quality  raw  ultrasound  data  as  a  basis  for  the  subsequent 
work  elements.  This  transducer  is  operable  at  up  to  14MHz  whereas  the  previously  available 
transducer  was  only  operable  up  to  8  MHz.  Consequently,  we  obtain  excellent,  and  significantly 
improved,  image  resolution. 

Aim  2.  Develop  and  test  a  tissue  elasticity  imaging  system. 

We  assembled  the  apparatus  to  enable  the  new  approach  to  transrectal  ultrasound  based  strain 
imaging  in  Year  1  and  Year  2.  We  have  also  made  several  custom  prostate  phantoms  using 
locally  developed  techniques  (Negron,  Viola  et  al.  2002).  Using  internally  made  phantoms,  we 
have  been  able  to  iterate  efficiently  the  design  and  also  to  fabricate  replacement  phantoms  at 
low  cost  in  a  timely  manner.  Phantoms  generally  deteriorate  over  time  due  to  dehydration 
through  the  membranes.  We  have  tested  the  tissue  elasticity  system  using  both  an  older  8  MHz 
transrectal  transducer  and  the  newer  14  MHz  transducer  connected  to  our  Siemens  Sequoia 
ultrasound  machine.  (Our  techniques  can  be  migrated  to  other  ultrasound  systems  if  resource 
and  contractual  issues  are  addressed.) 

We  spent  a  considerable  portion  of  Year  3  working  on  improved  algorithms  for  3D  elastography 
and  image  processing.  This  is  described  below: 
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Multi-resolution  hybrid  strain  estimator 

Time  delay  estimation  has  been  extensively  reported  in  the  literature  as  an  effective  method  to 
estimate  blood  flow  (Kasai  et  al.  1985)  and  tissue  motion  (Chen  et  al.  1992).  In  the  field  of 
elastography,  O’Donnell  et  al.  (1994))  have  used  the  zero-crossing  of  the  phase  of  the  complex 
cross  correlation  function  (base  band  signal)  to  estimate  tissue  strain,  whereas  Konofagou  et  al. 
(2000a)  have  used  shift  in  the  power  spectrum  of  the  pre-  and  post-compression  echo  RF  signal 
to  estimate  tissue  strain.  Bilgen  (1999)  has  reported  a  wavelet  based  direct  strain  estimator  (no 
gradient  operation)  for  elastography.  In  this  section,  we  describe  a  multistage  hybrid  adaptive 
displacement  (and  strain)  estimation  algorithm.  Figure  1  illustrates  the  complete  flow  chart  of 
this  strain  estimation  algorithm,  which  is  elaborated  in  the  following  paragraphs.  In  strain 
estimation,  cross-correlation  between  pre-  and  post-compression  signals  is  maximized  when 
two  signals  are  jointly  stationary.  In  addition,  if  the  time-bandwidth  product  is  large  and  the  shifts 
in  the  signal  are  small,  cross-correlation  is  an  efficient  estimator  of  displacement  (hence  strain) 
and  in  the  limit  achieves  the  Cramer-Rao  lower  bound  on  variance  (Walker  and  Trahey  1995). 
However,  in  practice,  pre-  and  post-compression  signals  are  jointly  non-stationary.  Hence, 
signal  companding  (Chaturvedi  at  al.  1999)  or  stretching  (Cespedes  and  Ophir  1993)  is 
frequently  used  to  improve  image  quality.  However,  for  global  stretching  to  be  effective,  an  a 
priori  estimate  of  applied  tissue  strain  is  required.  In  conventional  elastography,  a  calibrated  and 
computer  controlled  motion  stage  is  used  to  apply  an  external  displacement,  it  accurately 
determines  the  applied  displacement  at  the  top  surface  of  the  object  being  compressed  due  to 
uniform  compression.  Hence,  the  applied  strain  at  the  top  of  the  compressed  object  is  known  a 
priori.  In  hand-held  and  trans-rectal  prostate  elastography,  applied  displacement  cannot  be 
guessed  a  priori  as  the  applied  compressions  are  not  calibrated  or  controlled.  In  addition,  the 
strain  induced  in  the  compressed  object  decays  or  attenuates  with  depth,  and  is  governed  by 
non-  uniform1  boundary  conditions. 


1  When  the  compressor  plate  is  smaller  in  dimension  than  that  of  the  compressed  homogeneous  object,  the  stress  distribution  in  the  object  is 
non-uniform  and  is  referred  to  as  non-uniform  boundary  condition  at  the  object  boundary  (Shapo  et  al.  1996,). 
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Figure  1:  Flow  chart  of  the  new  strain  estimation  algorithm.  Note:  dotted  lined- 
rectangular  boxes  denote  first-stage  (I),  whereas,  dashed  lined-rectangular  boxes  signify 
second-stage  (II),  and  solid  lined-rectangular  boxes  represent  third-stage  (III). 

Hence,  uniform  global  stretching  in  the  first  stage  of  any  adaptive  strain  estimation  algorithm 
(Srinivasan  et  al.  2002b)  may  produce  sub-optimal  results.  Also,  when  applied  strain  is  large  (> 
2%,  for  a  ~4  cm  deep  prostate  tissue  imaged  at  12  MHz),  the  induced  displacement  is 
sufficiently  larger  than  the  wavelength  corresponding  to  center  frequency  and  may  introduce 
ambiguity  in  peak  detection  of  the  cross-correlation  function  obtained  by  tracking  the  pre-  and 
post-compression  RF  data.  Hence,  in  the  first  stage  of  our  algorithm,  we  demodulate  the  RF 
data  and  track  the  envelope  of  the  pre-  and  post-compression  echo  RF  frames  to  yield  average 
’’first-guess”  strain  estimates  over  the  tracking  window-W.  The  window,  W,  is  the  spatial  length 
of  the  kernel  over  which  the  pre-  and  the  post-compression  signal  cross-correlation  is  carried 
out.  The  window  in  this  stage  is  larger  than  in  the  subsequent  RF  tracking  stages  to  obtain 
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estimates  with  higher  SNR.  The  displacement  obtained  from  this  stage  is  referred  to  as  the  first 
stage  axial  displacement  (1Da),  where  superscript  ‘T  denotes  the  stage  and  subscript  ‘a’ 
denotes  the  displacement  type  (axial 

In  the  second  stage,  a  polynomial  curve2  is  fitted  to  the  estimated  local  strain  profile  over 
an  A-line,  and  the  corresponding  A-line  is  stretched  by  the  curve-fitted  first-stage  strain. 
Srinivasan  et  al.  (2002b)  stretched  the  post-compression  data  globally  by  the  average  applied 
strain,  whereas  Chaturvedi  et  al.  (1999)  compressed  the  pre-compression  data  by  the  average 
tissue  strain.  In  either  case,  the  first  stage  did  not  take  into  account  the  local  variations  in  tissue 
contrast.  In  addition,  both  algorithms  require  an  a  priori  estimate  of  the  applied  strain.  In  the 
proposed  algorithm,  signal  stretching  is  performed  locally  over  a  tracking  window  on  every  A- 
line,  taking  account  of  the  local  variation  in  tissue  elasticity  without  any  a  priori  knowledge  of  the 
applied  tissue  strain.  Lateral  displacement  (2D|)  is  estimated  by  performing  lateral  tracking 
(Konofagou  et  al.  2000b).  The  pre-  and  the  stretched  post-compression  echo  RF  data  are  then 
tracked  to  yield  a  local  strain  estimate  over  the  search  window.  The  axial  displacement 
estimated  in  this  stage  is  added  to  the  first  stage  axial  displacement  and  the  resultant  axial 
displacement  is  referred  to  as  the  second  stage  axial  displacement  (2Da). 

The  third  stage  is  used  to  estimate  residual  strain.  In  the  third  stage,  the  original  un-stretched 
post-compression  RF  data  are  stretched  locally  over  a  tracking  window  using  the  estimated 
second-stage  strain.  The  stretched  post-compression  data  are  then  tracked  with  the  pre¬ 
compression  data  to  yield  residual  axial  displacements.  The  lateral  displacements  are  obtained 
by  tracking  the  A-lines  laterally.  When  performing  lateral  tracking,  the  search  window  in  the  pre¬ 
compression  data  is  tracked  with  the  corresponding  window  on  adjacent  A-lines  in  the  post¬ 
compression  data.  The  maximum  cross-correlation  coefficients  over  these  windows  are 
estimated  and  a  cosine  curve  (cosine  interpolation)  (Cespedes  and  Ophir  1993)  is  fitted  to  the 
data  points.  The  shift  in  the  peak  of  this  curve  in  the  transverse  (lateral  or  azimuthal)  direction 
yields  the  lateral  displacement  of  the  local  post-compression  window.  The  axial  and  lateral 
displacements  estimated  in  this  stage  are  added  to  the  second  stage  axial  and  lateral 
displacement  to  yield  the  final  axial  (3Da),  and  lateral  displacements  (3D|).  The  gradient  of  the 
axial  displacement  image  yields  a  strain  image,  but  gradient  operator  is  a  high  pass  operator 
and  amplifies  noise.  A  least  square  estimator  may  be  applied  to  a  displacement  image  to 
produce  a  smooth  elastogram,  but  the  losses  in  tissue  contrast  and  resolution  may  be 
substantial  (Kallel  and  Ophir  1997).  A  staggered  strain  operator  (Srinivasan  et  al.  2002a)  limits 
noise  amplification  and  provides  smooth  elastograms  while  minimizing  the  loss  in  contrast  and 
resolution.  Hence,  in  the  last  step  of  this  algorithm,  we  apply  staggered  strain  estimation 
technique  to  the  final  stage  axial  displacement  image  to  generate  a  final  stage  elastogram.  In 
principle,  the  staggering  estimator  performs  a  band-pass  filtering  operation  partially  mitigating 
the  noise-enhancing  effect  of  the  high  pass  filtering  resulting  from  the  gradient  operator. 


2  The  estimated  first  stage  strain  profile  is  interpolated  using  cosine  interpolation  to  match  the  sampling  rate  of  the  RF  data. 


3D  Elastographic  simulation  framework 


An  elastographic  simulation  framework  similar  to  that  proposed  by  Nightingale  et  al.  (2000)  and 
Patil  et  al.  (2006)  was  used  to  in  this  work.  Patil  et  al.  had  used  a  3D  shift-invariant  convolution 
model  for  generating  the  pre-  and  post  compression  echo  RF  data.  In  this  work,  we  have  used  a 
more  realistic  Field  II  simulation  model  (Jensen  and  Nikolov  2000)  for  echo  RF  data  generation. 
Most  FEA  simulations  performed  by  Patil  et  al.  were  restricted  by  the  assumption  of  uniform 
boundary  conditions,  which  is  not  strictly  valid  for  transrectal  prostate  elastography.  In  this  work, 
all  FEA  simulations  were  performed  with  non-uniform  boundary  conditions.  Figure  2  illustrates 
the  flow  diagram  of  the  simulation  framework  used  in  this  work.  ANSYS  (Canonsburg,  PA),  a 
finite  element  analysis  (FEA)  simulation  software  package,  was  used  in  this  work.  Arbitrary 
regions  of  interest  (ROI’s)  can  be  created  within  the  FEA  model  by  changing  the  local  material 
properties  (shear  modulus,  material  density  etc.).  A  fine  FEA  mesh  (approximately  13  FEA 
nodes  per  resolution  cell)  results  in  accurate  parameter  estimation  (3D  internal  displacement,  in 
this  case).  Boundary  conditions  (applied  displacement)  can  be  applied  to  an  object  surface,  an 
object  element,  or  a  node  of  the  object  element  (after  meshing).  The  internal  object  deformation 
due  to  the  applied  displacement  is  then  estimated  by  solving  the  associated  partial  differential 
equations  numerically.  In  this  work,  a  preconditioned  conjugate  gradient  iterative  solver  was 
chosen  to  perform  numerical  simulations  (Elman  1981).  This  solver  was  chosen  because  it 
provided  similar  performance  as  compared  with  the  other  solvers  (within  ANSYS)  at  a  relatively 
small  computational  time  (Patil  2005).  A  grid  of  scatterers  was  defined  over  the  3D  ROI  of  the 
simulated  object  (48  scatterers  per  resolution  cell  (Walker  2006)).  The  post-compression 
scatterer  position  was  simulated  by  applying  the  3D  displacement  field  obtained  from  the  FEA 
software  to  the  scatterers. 


Figure  2:  3D  elastographic  simulation  framework  . 

It  is  worthwhile  noting  that  the  FEA  simulations  accounted  not  only  for  translational  but  also  for 
rotational  motion  of  an  object  due  to  any  arbitrary  boundary  condition.  Thus,  realistic  tissue 
motions  can  be  simulated  for  arbitrary  boundary  conditions  and  arbitrary  object  shapes.  Field  II 
simulations  were  used  to  generate  the  pre-  and  post-compression  echo  RF  data.  The  generated 
pre-  and  post-compression  data  can  be  tracked  using  different  displacement  (strain)  estimation 
algorithms  to  generate  a  strain  image.  In  this  work,  we  used  the  strain  estimation  (signal 
processing)  algorithm  detailed  in  section  2. 


Performance  evaluation  of  the  strain  estimator  using  3D  elastographic  simulations  and 
experiments 

A.  Simulations 


Elastographic  simulations  were  performed  using  the  framework  described  in  section  3.  Table  1 
lists  the  simulation  parameters  used  for  performing  various  simulations  described  in  this  section. 

Table  1  Simulation  parameters  . 
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Entity 


Value 


Dimensions  of  background 

Modulus  ratio  of  inclusion  and  background 

Approximate  diameter  of  the  inclusion 

Center  frequency  of  scan 

Sampling  frequency 

System  SNR 

ID  linear  array 

Focus 

Fractional  bandwidth 
Window  length  (W) 

Overlap 


40x40x40  mm  cube 
3 

10  mm 
12  MHz 
80  MHz 
30  dB 

192  elements,  0.2  mm 

pitch 

20  mm 

40% 

2  mm 
90% 


A  cube  (40  x  40  x  40  mm)  enclosing  a  spiculated  inclusion  (10  mm  diameter)  was  used  as  a 
mechanical  model  to  simulate  the  background  and  the  lesion  (Figure  3).  The  embedded 
inclusion  (Young’s  Modulus=12  KPa)  was  simulated  to  be  three  times  stiffer  than  the 
background  (Young’s  Modulus=4  KPa).  A  uniform  cube  was  used  for  the  SNR  study.  The  cubes 
were  compressed  axially  from  the  top  and  assumed  to  be  rested  on  a  fixed  surface  such  that 
their  motions  in  the  lateral  and  elevational  directions  were  unconstrained  (slip  conditions).  The 
cubes  were  compressed  along  the  elevational  axis  of  symmetry  using  a  compressor  plate  (40  x 
3  mm-  lateral  elevation  plane)  and  the  transducer  array  was  placed  over  the  compressor  plate. 
The  surface  area  of  the  compressor  plate  was  smaller  than  the  surface  area  of  the  top  of  the 
cube  to  simulate  non-uniform  boundary  conditions  at  the  top  surface  of  the  cube  (Figure  3).  This 
results  in  non-uniform  stress  distribution  through  the  depth  of  the  tissue  and  consequently 
causes  strain  decay.  The  cube  was  subjected  to  a  range  of  applied  strain  (0.5-10%).  ANSYS 
was  used  to  perform  these  simulations.  The  transducer  was  focused  at  20  mm  (range  or  axial 
direction).  RF  data  were  generated  using  the  3D  elastographic  simulation  framework  as 
described  in  section  3 


Figure  3:  Schematic  of  the  simulated  phantom  and  imaging  array  used  for  CNRe  study. 
The  cube  is  40x40x40  mm  and  the  compressor  plate  is  40x  3  mm. 
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The  generated  pre-  and  post-compression  RF  data  were  tracked  using:  1)  the  proposed 
algorithm,  2)  the  algorithm  proposed  by  Srinivasan  et  al.  (2002b),  and  c)  ID  RF  data  tracking 
using  uniform  global  stretching  (Cespedes  and  Ophir  1993).  Signal-to-noise  ratio  (O’Donnell  et 
al.  1994)  and  contrast-to-noise  ratio  (Bilgen  and  Insana  1997)  were  used  as  metrics  for 
comparing  the  above  listed  tracking  algorithms.  Figures  4  through  7  illustrate  the  first,  second 
and  final  stage  elastograms  obtained  using  the  proposed  algorithm.  Figure  8  illustrates  the 
numerical  strain  image  obtained  from  3D  FEA  simulation 


Figure  4:  First  stage  elastogram  (Left)  and  axial  displacement  image 
(Right).  The  field  of  view  is  40x40  mm.  The  images  were  obtained 
using  parameters  specified  in  Table  1. 


Figure  5:  Strain  profile  over  an 
A-line  in  the  elastogram  (Figure 
4-Left)  obtained  by  the  B-mode 
tracking. 


Figure  6:  Second  stage  elastogram  (Left),  axial  displacement  image  (Center),  lateral  displacement 
image  (Right).  The  field  of  view  is  40x40  mm.  The  images  were  obtained  using  parameters  specified  in 
Table  1. 


Figure  7:  Final  stage  elastogram  (Left),  axial  displacement  image  (Center),  lateral  displacement  image 
(Right).  The  field  of  view  is  40x40  mm.  The  images  were  obtained  using  parameters  specified  in  Table 
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Figure  8:  Ideal  strain  image  obtained  from 
the  FEA  simulations.  The  image  was 
obtained  using  non-uniform  boundary 
conditions,  other  mechanical  parameters 
used  in  FEA  simulations  are  listed  in  Table 
1. 
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Figure  9  illustrates  the  theoretical  elastographic  SNR  (SNRe)  plot  obtained  using  the 
specifications  of  the  simulated  transducer  and  the  aforementioned  signal  processing 
parameters.  The  theoretical  framework  proposed  by  Varghese  and  Ophir  (1997)  was  used  for 
generating  the  theoretical  SNRe  curve  plot.  It  is  important  to  note  that  the  theoretical  SNRe 
curve:  1)  assumes  one  dimensional  speckle  motion,  and  2)  joint  stationarity  of  pre-  and  post¬ 
compression  signals,  3)  neglects  frequency  dependent  attenuation  and  4)  assumes  that  the 
observation  window  is  at  the  focus  of  the  transducer.  In  view  of  these  assumptions,  the 
theoretical  SNRe  can  never  be  achieved  in  practice.  Nevertheless,  the  theoretical  SNRe  plot  can 
be  used  as  a  “gold  standard”  to  evaluate  the  performance  of  any  strain  estimation  algorithm. 
Figure  10  illustrates  the  signal-to-noise  curves  estimated  using  the  proposed  algorithm  and  the 
above  listed  tracking  approaches.  Equation  (1)  was  used  to  estimate  the  elastographic  signal- 
to-noise  ratio  (SNRe)  from  the  reconstructed  elastograms  (O’Donnell  et  al.  1994), 


where  s  is  the  average  strain  estimated  over  the  entire  image.  os  is  the  variance  of  strain 
estimates.  The  SNRe  curve  estimated  using  the  proposed  algorithm  envelopes  the  SNRe  curves 
obtained  using  the  global  stretching-based  RF  data  tracking  algorithm  and  the  algorithm  by 
Srinivasan  et  al.  (2002b).  For  all  values  of  applied  strain,  the  SNRe  plot  obtained  from  the 
proposed  algorithm  was  statistically  different  (p  <  0.05)  from  the  one  obtained  using  the 
algorithm  proposed  by  Srinivasan  et  al.  (2002b).  Analysis  of  variance  was  used  as  a  statistical 
measure  for  quantifying  the  differences  between  the  two  SNRe  curves.  Thirty  independent 
realizations  were  used  for  the  statistical  study.  Figure  11  illustrates  the  image  variance 
estimated  at  different  values  of  applied  strain.  As  predicted  using  strain  filter  theory  (Varghese 
and  Ophir  1997),  the  variances  in  strain  estimation  increase  non-linearly  with  applied  strain.  The 
rate  of  increase  in  variance  is  lowest  when  using  the  proposed  algorithm.  Figure  12  illustrates 
the  average  cross-correlation  coefficient  (p)  over  the  image  obtained  using  different  tracking 
algorithms.  The  cross-correlation  coefficient  obtained  by  tracking  the  pre-  and  the  post¬ 
compression  signal  can  be  expressed  analytically  as: 


(2) 


(3) 


(4) 

(5) 
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1  I  A 

<j2P=  | S2  P{t)dt 

- T 12  ,  and 

772 

o'2  s  =  \s2s(t)dt. 

-77 2 

Subscript  p  denotes  pre-compression  and  s  denotes  post-compression  (O’Donnell  et  al.  1994). 
T  is  the  tracking  window  (W=  T*c/2)\  c  is  the  speed  of  sound  through  the  tissue.  pps  is  the  cross¬ 
correlation  coefficient.  The  error  bars  in  Figure  12  illustrate  the  variance  in  the  average  image 
cross-correlation  coefficient  obtained  from  thirty  independent  realizations.  The  rate  of  decrease 
in  the  correlation  coefficient  as  a  function  of  applied  strain  is  smallest  for  the  proposed 
algorithm.  Figure  13  illustrates  different  contrast-to-noise  ratio  curves-  CNRe  (Bilgen  and  Insana 
1997)  obtained  from  the  elastograms  simulated  using  different  strain  estimation  algorithms. 
Contrast-to-noise  ratio  is  expressed  as, 

cm.- ^4  <6) 

I*2  t+(Jb2 

where,  st  is  the  average  strain  in  the  target,  and  Sb  is  the  average  strain  of  the  background.  o2tis 
the  variance  of  the  strain  estimates  in  the  target,  and  o2bis  the  variance  of  the  strain  estimates 
in  the  background.  CNRe  values  for  a  range  of  applied  strain  (0.5-10%)  obtained  using  the 
proposed  algorithm  are  3  dB  above  the  CNRe  curve  obtained  using  the  algorithm  proposed  by 
Srinivasan  et  al. (2002b)  and  are  6  dB  higher  than  the  RF  data  tracking  algorithm  with  uniform 
global  stretching  (Cespedes  and  Ophir  1993). 


10'  10'  10  10  10 
Applied  Strain  (%) 

Figure  9:  Theoretical  SNRe  plot  for  a  range  of  applied  strains.  The  plot  was  obtained  at  12 
MHz  center  frequency,  40  %  fractional  bandwidth  (-6  dB),  system  SNR  of  30  dB,  window 
length  of  2  mm  and  overlap  of  90  %. 
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Figure  10:  SNRe  plots  comparing  proposed  Figure  11:  Estimated  image  variance 
algorithm  with  uniform  global  stretching  comparing  proposed  algorithm  with 
only  (Cespedes  and  Ophir  1993),  and  global  uniform  global  stretching  only  (Cespedes 
followed  by  local  stretching  approaches  and  Ophir  1993),  and  global  followed  by 
(Srinivasan  et  al.  2002b).  local  stretching  approaches  (Srinivasan  et 


al.  2002b). 
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Figure  12:  Average  cross-correlation 
coefficient  over  the  reconstructed  image 
obtained  using  proposed  algorithm, 
algorithm  with  uniform  global  stretching 
only  (Cespedes  and  Ophir  1993),  and  global 
followed  by  local  stretching  approaches 
(Srinivasan  et  al.  2002b). 
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Figure  13:  CNRe  ratio  plots  obtained  using 
proposed  algorithm,  algorithm  with 
uniform  global  stretching  only  (Cespedes 
and  Ophir  1993),  and  global  followed  by 
local  stretching  approaches  (Srinivasan  et 
al.  2002b). 
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B.  Experiments 


Experimental  evaluation  of  the  tracking  algorithm  was  performed  on  an  acryl-amide  based 
prostate  phantom  previously  fabricated  by  (Li  et  al.  2005).  The  mean  cross-sectional  diameter  of 
the  embedded  oblong  inclusion  in  the  longitudinal  (range)  direction  was  approximately  7  mm.  A 
Siemens  Sequoia  512  scanner  (Siemens  Medical  Solutions,  Mountain  View,  CA)  was  used  in 
this  study.  Multiple  demodulated  “l/Q”  (In-phase/Quadrature)  radio  frequency  beamformed  lines 
of  acoustic  data  were  acquired  from  the  ultrasound  scanner  using  a  research  interface 
employing  an  l/Q  data  capture  board.  A  192  element,  12  MHz,  -6  dB  40  %  fractional-bandwidth 
transrectal  ultrasound  transducer  was  used  for  scanning  the  prostate  phantom.  The  depth  of 
acquisition  was  24  mm  and  the  width  of  acquisition  was  40  mm.  The  experimental  apparatus 
and  the  protocols  used  for  the  validation  are  described  in  (Li  et  al.  2005)  and  will  be  reiterated  in 
detail  in  the  phantom  experiment  section  (Section  6).  Fifty  frames  were  acquired  such  that  the 
total  applied  strain  over  50  frames  was  2  %.  Every  25th  frame  (1  and  25-E125,  2  and  26-E2,26 
and  so  on...)  was  tracked.  Final  elastogram  was  generated  by  averaging  over  five  independently 

obtained  elastograms  ((E125  + . +  E5  29)/5)  Figure  14a  illustrates  an  elastogram  obtained  using 

the  method  proposed  in  this  paper.  Figure  14b  illustrates  an  elastogram  reconstructed  by  using 
the  algorithm  proposed  by  Srinivasan  et  al.  in  (2002b).  As  illustrated  in  Figures  14a  and  14b, 
elastograms  generated  by  the  algorithm  proposed  in  this  paper  provide  higher  CNReand  retain 
the  shape  of  the  enclosed  inclusion.  The  depth  dependent  strain  attenuation  is  clearly  visible  in 
the  elastograms  in  Figures  14a  and  14b.  Thus  uniform  global  stretching  over  the  entire  image 
as  suggested  in  the  first  stage  of  the  algorithm  by  Srinivasan  et  al.  (2002b)  results  in  stretching 
inaccuracies  and  is  a  likely  cause  of  the  noise  observed  in  Figure  14b.  The  algorithm  proposed 
in  this  paper  stretches  every  A-line  by  the  local  strain  profile  over  that  A-line  as  predicted  by  the 
B-mode  tracking  stage  (the  first  stage  of  the  proposed  algorithm).  Also,  the  stretching  takes  into 
account  the  strain  decay  with  depth  and  hence  minimizes  noise  due  to  stretching  inaccuracies, 
which  are  likely  to  occur  due  to  tissue  inhomogeneity.  The  algorithm  by  Srinivasan  et  al.  (2002b) 
requires  an  a  priori  estimate  of  the  applied  strain  for  stretching  the  post-compression  data, 
whereas  the  algorithm  proposed  in  this  paper  does  not  have  such  a  requirement.  Thus,  for 
trans-rectal,  and  handheld  elastography,  the  proposed  algorithm  is  expected  to  exhibit  superior 
performance  to  algorithm  proposed  by  Srinivasan  et  al.  (2002b).  For  the  case  with  uniform 
boundary  condition  (as  in  motion  stage  based  elastography,  where  the  applied  strain  is  known  a 
priori),  we  expect  that  both  algorithms  should  demonstrate  similar  performance. 
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Figure  14a  The  modulus  contrast  ratio  of  the 
inclusion  in  the  image  with  respect  to  the 
background  was  10.  A  window  length  of  2 


Figure  14b:  The  modulus  contrast  ratio  of 
the  inclusion  in  the  image  with  respect  to 
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mm  and  an  overlap  of  90  %  were  used  for  the  background  was  10.  A  window  length  of 
processing  the  elastogram.  The  obtained  2  mm  and  an  overlap  of  90  %  were  used  for 
elastogram  is  an  average  strain  image  over  processing  the  elastogram.  The  obtained 
5  frames.  The  elastogram  was  reconstructed  elastogram  is  an  average  strain  image  over 
using  the  proposed  algorithm.  5  frames.  The  elastogram  was  by 

reconstructed  using  the  algorithm  proposed 
by  Srinivasan  et  al.  (2002b). 
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1.  Rendering  3D  shapes  of  embedded  inclusions  :Simulations 

3D  elastographic  simulations  framework  described  in  section  3  was  used  for  3D  shape 
rendering  study.  Two  different  cases  were  considered  for  the  3D  rendering  study-  the  axial  and 
the  lateral  resolution  study  (Figure  15).  The  transducer  array  (Field  II)  was  swept  across  (along 
elevation)  the  inclusions  to  generate  elastograms  by  3D  rendering.  The  spiculated  inclusion 
models  a  malignant  lesion,  whereas  the  smooth  spherical  inclusion  models  a  benign  lesion 
(Ueno  1986,  Lee  1988).  In  the  model  shown  in  Figure  15,  the  background  was  simulated  as  a 
40x40x40  mm  cube  with  uniform  density  and  stiffness  of  4  KPa  The  inclusion  was  three  times 
as  stiff  as  the  surrounding  background.  The  spherical  inclusion  had  a  cross-sectional  diameter 
of  5  mm;  whereas  the  spiculated  malignant  inclusion  had  a  mean  cross-sectional  diameter  of 
approximately  6  mm.  External  displacement  was  applied  to  the  top  surface  of  the  cube  as 
illustrated  in  Figure  15.  Parameters  listed  in  Table  1  were  used  to  for  simulations  described  in 
this  section.  The  simulated  transducer  was  slowly  swept  across  the  top  surface  of  the  cube  and 
twenty  three  pre-  and  post-compression  RF  data  slices  were  collected  along  the  elevational 
dimension.  The  acquired  pre-  and  post-compression  RF  data  sets  were  tracked  using  the 
algorithm  detailed  in  section  2.  Various  elastograms  obtained  along  the  elevational  direction 
were  segmented  using  3D  gradient  vector  flow  active  surface  (GVF-AS)  algorithm  (Tay  et  al. 
2006)  and  rendered  in  3D  to  reconstruct  the  3D  shapes  of  the  inclusions.  The  strain  images 
obtained  directly  from  the  FEA  simulations  were  also  segmented  and  rendered  to  produce  ideal 
3D  shapes  of  the  modeled  inclusions.  Figure  16  illustrates  the  3D  reconstruction  from  the  axial 
resolution  study.  The  modeled  malignant  and  benign  lesions  were  embedded  in  the  3D  cube 
and  were  placed  longitudinally  along  the  axis  of  symmetry.  The  axial  separation  between  the 
two  lesions  was  reduced  in  incremental  steps  of  0.1  mm.  For  every  step,  the  lesions  were 
reconstructed  from  the  elastograms  and  rendered  in  3D.  The  process  was  repeated  until  the 
lesions  were  close  enough  to  each  other  such  that  their  reconstruction  resulted  in  one  fused 
entity  instead  of  two  distinct  entities.  The  axial  resolution  of  the  3D  reconstruction  was 
approximately  0.8  mm  (five  wavelengths).  Figure  16-a  illustrates  a  3D  reconstruction  from  the 
generated  elastograms,  whereas  Figure  16-b  illustrates  a  3D  reconstruction  from  the  ideal 
numerical  strain  images  obtained  from  the  FEA  simulations.  As  expected,  the  3D  reconstruction 
from  ideal  images  is  sharper  than  that  from  the  elastogram.  Also,  the  lesions  reconstructed  from 
the  elastogram  exhibit  shape  warping  that  can  be  attributed  to  limitations  arising  from  the  signal 
processing  parameters  (window  length  and  window  overlap  etc.),  and  the  non-uniformity  of  the 
sampling  of  the  ultrasound  data  in  the  three  dimensions  (Patil  et  al.  2006).  Lateral  resolution  of 
the  3D  reconstruction  was  estimated  by  repeating  the  steps  detailed  for  the  axial  resolution 
study.  The  lateral  resolution  of  the  3D  reconstruction  was  approximately  1.5  mm  (ten 
wavelengths).  Figures  17  a-b  illustrate  a  3D  reconstruction  from  the  lateral  resolution  study 
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Figure  15:  The  a  ive  figures  illustrate  the  two  different  Simula  b  i  cases  considered  for 
the  3D  rendering  study.  The  arbitrary  shaped  inclusion  moaeis  a  malignant  lesion, 
whereas  the  smooth  spherical  inclusion  models  a  benign  lesion  (Ueno  1986,  Lee  1988). 
The  cases  a  and  b  were  considered  to  evaluate  the  axial  and  the  lateral  resolution  of  the 
3D  reconstruction. 
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Figure  16:  Axial  resolution  simulation  study:  (a)  3D  reconstruction  from  the 
reconstructed  elastograms,  (b)  3D  reconstruction  from  the  numerical  strain  images 
derived  from  the  FEA  simulations 
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Figure  17:  Lateral  resolution  simulation  study:  (a)  3D  reconstruction  from  the 
reconstructed  elastograms,  (b)  3D  reconstruction  from  numerical  strain  images 
derived  from  the  FEA  simulations. 
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2.  Elastographic  phantom  experiments  for  3D  rendering  study 

A.  A  purpose  built  prostate  phantom 

The  acrylamide-based  prostate  phantom  was  fabricated  using  the  protocol  published  by  Negron 
et  al.  (2002).  The  prostate  phantom  consisted  of  three  compartments  -  the  prostate  tissue  with 
embedded  inclusion,  the  anal  passage,  and  the  surrounding  viscera.  A  schematic  of  the 
constructed  phantom  is  illustrated  in  Figure  18,  whereas  Figure  19  illustrates  a  photographic 
view  of  the  prostate  phantom. 


Figure  18:  Schematic  representation  of  the  Figure  19:  Prostate  phantom-  a 
prostate  phantom.  1:  Anal  passage,  2:  photographic  view. 

Prostate  tissue,  3:Viscera 


The  prostate  material  and  the  surrounding  gel  (viscera)  were  composed  of  a  4%  acrylamide 
solution,  while  inclusions  in  the  prostates  were  composed  of  8%  acrylamide,  which  imparted 
different  stiffness  to  the  lesion  and  the  surrounding  prostate  tissue.  The  Young’s  modulus  for  4 
%  gel  is  approximately  4  KPa,  whereas  it  is  16  KPa  (Konofagou  et  al.  2003)  in  the  8%  gel  stiff 
inclusions.  Thus,  the  modulus  ratio  for  all  embedded  lesions  was  approximately  four.  Sephadex 
(GE-Amersham,  Piscatway,  NJ)  was  added  to  the  prostates  and  inclusions  to  give  speckle;  a 
higher  concentration  of  sephadex  was  added  to  the  inclusions  than  the  surrounding  prostates, 
yielding  a  different  speckle  pattern.  Two  prostates  of  approximately  100  ml_  each  were  included 
in  the  phantom,  one  containing  a  spherical  inclusions  and  the  other  containing  an  irregular 
inclusions.  The  smooth  inclusion  was  approximately  8  mm  in  diameter  whereas  the  irregular 
shaped  inclusion  had  a  mean  diameter  of  1  cm.  To  ensure  uniform  speckle-pattern  in  inclusions 
and  prostate,  the  mixture  of  different  component  solutions  was  enclosed  in  a  latex  balloon  and 
was  rolled  gently  for  a  setting  period  (approximately  20  minutes  following  the  mixing  of 
components)  until  the  solution  had  gelled.  The  rolling  protocol  mitigates  against  the  settling  of 
the  embedded  inclusion  and  may  sometimes  lead  to  a  situation  where  two  inclusions  are  placed 
apposing  each  other,  hence  two  different  prostates  were  constructed,  each  embedding  one 
inclusion. 

B.  “  I-Beam”  transducer  and  elastographic  phantom  experiments 

A  Siemens  Sequoia  512  scanner  was  used  for  image  acquisition.  The  prostate  phantom  was 
scanned  using  an  "I-Beam"  type  transducer  -  a  modified  linear  transducer  array  consisting  of  a 
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central  “imaging”  array  bordered  by  two  perpendicular  “tracking”  arrays  in  a  cylindrical  housing, 
as  illustrated  in  Figures  20  a-b.  Using  this  configuration,  the  angular  separation  between  two 
successive  image  slices  could  be  tracked  by  measuring  the  azimuthal  motion  in  the  tracking 
arrays  (Hossack  et  al.  2002,  Garson  et  al.  2006).  The  arrays  were  operated  at  12  MHz,  using  a  - 
6  dB  fractional  bandwidth  of  approximately  40  %.  A  small  tube  for  water  delivery  was  attached 
to  the  transducer  housing,  and  the  assembly  was  surrounded  with  a  water-tight  latex  transducer 
cover.  When  acquiring  different  image  slices,  the  transducer  had  to  be  swept  at  an  angle  across 
the  region  of  interest  (ROI).  This  assembly  was  mounted  so  that  during  the  angular  rotation, 
translational  motion  was  restricted,  facilitating  pure  rotational  motion  of  the  transducer.  The 
mounted  transducer  was  inserted  into  the  prostate  phantom  and  used  to  perform  a  series  of 
two-dimensional  elastographic  scans.  At  every  angular  position,  the  phantom  was  compressed 
by  injecting  water  in  the  transducer  sheath  via  the  connected  tube  using  a  syringe  pump  while 
simultaneously  collecting  fifty  l/Q  data  frames.  For  each  lesion,  gain,  field  of  view,  and  focal 
depth  were  optimized  for  maximum  visibility  of  that  particular  lesion.  The  data  were  processed 
offline  as  detailed  in  section  2  to  generate  elastographic  image  slices. 


Cable 
lo  System 


Figure  20:  a)  Configuration  of  the  modified  transducer  comprising  a  central  imaging  array 
bordered  by  two  perpendicular  tracking  arrays.  During  a  scan,  the  transducer  is  rotated 
about  its  central  axis,  as  indicated,  b)  A  photographic  view  of  prostate  probe  with 
tracking  and  imaging  arrays. 

C.  Experimental  results 

Elastographic  experiments  were  conducted  as  described  in  sub-section  B.  Both  irregular 
shaped  and  smooth  inclusions  were  imaged.  Sixteen  elastographic  and  sonographic  image 
slices  were  acquired  by  sweeping  the  trans-rectal  prostate  transducer  across  the  spatial  extent 
of  both  inclusions.  Figure  21  illustrates  elastograms  and  corresponding  sonograms  for  a  few 
angular  positions  of  the  trans-rectal  prostate  transducer  while  sweeping  across  the  irregular 
shaped  inclusion.  Figure  22  illustrates  similar  images  for  the  smooth  shaped  inclusion.  Figure 
23  illustrates  images  of  the  smooth  lesion  reconstructed  from  sonograms  and  elastograms, 
respectively.  The  shape  and  the  size  of  the  smooth  inclusion  reconstructed  from  the  sonograms 
and  elastograms  are  approximately  similar;  the  smooth  inclusion  reconstructed  from 
elastograms  tapers  at  the  longitudinal  edges,  and  is  probably  due  to  the  low 
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Figure  21:  Few  slices  across  the  embedded  inclusion  mimicking  a  malignant  lesion.  Top 
row  illustrates  elastograms,  bottom  row  illustrates  sonograms.  The  field  of  view  is  30  x  30 
mm. 


Figure  22:  Few  slices  across  the  embedded  inclusion  mimicking  a  benign  lesion.  Top 
row  illustrates  sonograms,  bottom  row  illustrates  elastograms.  The  field  of  view  is  40  x 
30  mm. 
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Figure  23:  a)  Smooth  inclusion  mimicking  a  benign  lesion  reconstructed  from 
elastograms.  b)  Smooth  inclusion  reconstructed  from  sonograms. 
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pass  filtering  or  smoothing  effect  of  the  processing  windows  used  for  generating  elastograms 
(Srinivasan  et  al  2002a,).  Figure  24  illustrates  the  irregular  shaped  lesion  reconstructed  from  the 
elastograms  and  sonograms,  respectively.  The  shape  of  the  irregular  shaped  inclusion,  when 
reconstructed  from  the  elastograms,  has  blunt  edges  and  appears  larger  in  size  than  that 
reconstructed  from  the  sonogram  or  B-mode  image.  Similar  results  have  been  reported  by 
Garra  et  al.  (1997)  where  the  authors  demonstrate  the  same  by  comparing  the  sonographic  and 
elastographic  images  for  an  irregular  shaped  malignant  cancer  in  in-vivo  human  scans.  Thus, 
we  contend  that  the  3D  size  and  the  shape  of  the  reconstructed  lesion  (inclusion)  have  the 
potential  to  assist  in  discriminating  benign  cancers  from  malignant  cancers. 
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Figure  24:  a)  Irregular  shaped  inclusion  mimicking  a  malignant  lesion  reconstructed  from 
elastograms.  b)  Irregular  shaped  lesion  inclusion  reconstructed  from  sonograms. 
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Discussion 

In  elastography,  image  quality  is  contingent  upon  accurate  estimation  of  internal  tissue 
displacements.  Signal  companding  is  a  simple,  yet  elegant,  approach  of  improving 
elastographic  image  quality.  In  the  second  section  of  this  article,  we  proposed  an  adaptive 
stretching  based  tracking  algorithm.  The  algorithm  takes  into  account  the  local  variations  in 
tissue  elasticity  and  performs  iterative  local  stretching,  which  provides  a  superior  CNRe  to 
previously  proposed  adaptive  tracking  algorithms  (Srinivasan  et  al.  2002b).  When  encountering 
non-uniform  boundary  conditions,  as  in  prostate  elastography,  the  proposed  algorithm  is 
expected  to  perform  better  than  conventional  stretching  based  algorithms  as  they  require  an  a 
priori  estimate  of  the  applied  strain.  When  using  an  ultrasound  scanner  operating  at  higher 
frame  rates  with  a  high  bandwidth  transducer  and  performing  handheld  or  trans-rectal  prostate 
elastography,  the  net  motion  between  the  consecutive  image  frames  may  be  small.  Thus,  the 
correlation-based  tracking  in  the  first  two  stages  of  this  algorithm  may  be  replaced  by  a  sum-of- 
absolute-difference  approach  (Chaturvedi  and  Insana  1999),  whereas  the  third  stage  may  be 
replaced  by  temporal  tracking  of  zero  crossings  (Srinivasan  et  al. 2002b)  of  the  pre-  and  post¬ 
compression  echo  RF  signals.  Thus,  a  reduction  in  computation  time  may  be  achieved  by  an 
order  of  magnitude  without  any  substantial  degradation  in  image  quality.  Such  a  modified 
algorithm,  when  implemented  on  FPGA’s,  can  be  easily  incorporated  into  programmable 
scanners  such  as  the  ULTRASONIX  RP  (Richmond,  BC)  and  may  help  in  real  time  tissue 
elasticity  imaging.  In  a  PC  based  environment,  real-time  image  registration  and  rendering 
algorithms  (Tay  et  al.  2006)  may  be  used  to  render  multiple  image  slices  to  accurately 
reconstruct  a  3D  volumetric  elastogram.  For  this  work,  the  elevational  resolution  is  expected  to 
be  a  function  of  the  greater  of  the  two  properties  of  the  data  set,  the  elevational  beam  width  and 
the  angular  separation  between  two  consecutive  slices  and  intuitively  may  be  assumed  to  be  a 
linear  function  of  depth  as  the  separation  between  two  consecutive  slices  increases  linearly  as  a 
function  of  depth.  Thus,  for  detecting  and  rendering  smaller  lesions  deep  in  tissue,  it  may  be 
essential  to  have  sufficiently  high  angular  sampling  in  addition  to  having  a  higher  image  frame 
rate.  Nevertheless,  such  3D  elastography  imaging  systems  may  be  used  for  a  variety  of 
applications  such  as  real-time  volumetric  monitoring  of  HIFU  lesions  (Souchon  et  al.  2005),  real¬ 
time  volumetric  estimation  of  viscous  properties  of  tissue  using  non-linear  or  sinusoidal 
compression  techniques  between  consecutive  image  frames  thereby  imaging  the  temporal 
behavior  of  the  tissue  strain  (Erpelding  et  al.  2005),  and  temporal  monitoring  of  the  growth  of 
cancers. 
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Aim  3.  Apply  image  processing  algorithms  to  acquired  B-Mode  and  elasticity  images  to  improve 
the  quantification  of  detected  elastic  inhomogeneities 

Image  Processing  and  Quantification 

We  have  progressed  to  a  accurate  3D  surface  rendering  from  2D  slices  by  implementing  a  3D  gradient 
vector  flow  (GVF)  snake  algorithm  (Xu  and  Prince  1998).  The  2D  and  3D  GVF  snake  algorithm  relies  on 
the  edges  to  be  well  delineated  and  contrast  between  the  various  regions  to  be  well  defined  in  each  2D 
ultrasound  slice  to  produce  a  surface  that  resembles  the  actual  scanned  object.  In  ultrasound,  images 
are  affected  by  a  granular  interference  pattern  known  as  “speckle”.  Before  an  accurate  3D  surface 
rendering  can  be  attained, a  preprocessing  despeckling  step  is  needed  to  reduce  the  variances  in  pixel 
values  within  homogeneous  regions  while  contrast  between  distinct  regions  are  concurrently  enhanced. 
We  have  evaluated  a  wide  variety  of  well  known  methods  such  as  the  Nagao  and  Matsuyama  filter 
(Nagao  and  Matsuyama  1979),  the  Lee  filter(Lee  1980),  the  Frost  et  al.  filter  (Frost,  Stiles  et  al.  1982), 
the  Kuan  et  al.  filter  (Kuan,  Sawchuk  et  al.  1985),  the  adaptive  weighted  median  filter  proposed  by 
Loupas  et  a/.( Loupas,  McDicken  et  al.  1989),  the  Wiener  filter  (Kailath  1976)  SRAD  proposed  by  Yu  and 
Acton  (Yu  and  Acton  2002),  and  a  novel  stochastically  driven  method  design  specifically  for  3D  surface 
rendering  from  2D  slices  of  the  prostate  and  other  organs.  The  method  we  develop  specifically  for  the 
task  of  this  grant  is  a  stochastically  driven  compression  filter  called  the  “squeeze  box  filter”  (SBF).  Our 
quantitative  evaluation  using  a  Field  II  (Jensen  and  Svendsen  1992)  simulated  B  mode  ultrasound  image 
with  contrast  enhancement  performance  determined  by  a  modified  Fisher  discriminant  has  determined 
that  the  newly  developed  SBF  method  outperformed  the  other  methods  and  is  exceptional  in  providing 
the  needed  intra-region  reduction  in  variance  and  inter-region  contrast  enhancement  with  computational 
efficiency.  In  Figs.  25,  26,  and  27,  we  show  the  3D  surface,  the  side  view,  and  the  bottom  view, 
respectively,  of  the  rendering  we  attained  from  a  sequence  of  scans  we  acquired  of  an  egg  phantom. 

The  sequence  consists  of  acquiring  a  2D  slice  every  millimeter  along  the  long  axis  of  an  egg  phantom. 

We  processed  each  slice  with  the  SBF  despeckling  method.  The  3D  rendering  was  attain  by  SBF 
processing  each  slice  then  applying  a  3D  GVF  snake  to  attain  the  final  results  shown  in  Fig.  13,  14,  and 
15.  It  is  very  evident  that  our  method  captured  the  essential  size  and  shape  of  the  egg  phantom. 
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Figure  25.  3D  surface  found  by  the  3D  GVF  with  slices  preprocessed  by  SBF. 
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Figure  26.  Side  view  of  the  surface  rendered  by  the  3D  GVF  snake 
with  slices  preprocessed  by  SBF. 


Figure  27.  Bottom  view  of  the  surface  of  the  egg  phantom 
rendered  by  the  3D  GVF  snake  with  slices  preprocessed  by  SBF. 
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After  the  3-D  surface  of  the  prostate  is  segmented,  the  volume  is  determined  by  applying  a  novel 
blobbing  technique,  called  the  multi-directional  connected  component  analysis  (MDCCA).  The  MDCCA 
requires  that  the  surface,  which  is  represented  by  a  set  of  3-D  vertices  and  faces,  to  be  densely  meshed 
so  that  each  vertex  is  no  more  than  one  voxel  from  its  neighboring  vertices.  A  2-D  binary  image  of  the 
contour  is  defined  by  the  intersection  of  a  plane  through  the  densely  mesh  surface.  This  binary  image  is 
zero  except  when  the  pixel  lies  on  the  contour,  in  which  case  the  pixel  is  set  at  one.  The  interior 
enclosed  by  the  2-D  close  contour  is  determined  by  extracting  the  connected  component  in  several 
directions  and  applying  a  logical  “and”  operator  to  the  resulting  binary  images.  A  connected  component 
binary  image  in  some  fixed  direction  is  the  result  of  approaching  the  2-D  closed  contour  in  some  fixed 
direction,  say  left  to  right  in  each  row  or  top  to  bottom  in  each  column  for  examples,  and  setting  the  pixel 
values  to  zero  until  a  contour  pixel  is  reached,  in  which  case  the  remaining  pixels  in  the  fixed  direction 
are  set  to  one.  Fig.  28  shows  the  connected  component  binary  images  in  eight  directions.  Fig.  29 
shows  the  result  of  applying  a  logical  “and”  operator  to  the  eight  binary  images  in  Fig.  28.  The  MDCCA 
of  non-intersecting  slices  through  the  3-D  densely  meshed  surface  yields  a  binary  3-D  data  where  the 
voxels  are  one  when  the  voxel  lies  within  the  interior  of  the  3-D  surface  and  zero  otherwise.  An  example 
of  a  binary  3-D  data  produced  by  applying  the  MDCCA  to  non-intersecting  slices  is  shown  in  Fig.  30. 
The  volume  in  units  of  voxels  of  the  object  enclosed  by  the  3-D  surface  is  attained  by  summing  the 
binary  3-D  data  of  MDCCA  of  non-intersecting  slices.  The  volume  is  converted  to  units  of  cubic 
centimeter  (or  millimeter  or  other  units)  by  multiplying  with  the  voxel  resolution. 


Figure  29.  Connected  components  in  eight  directions. 
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Figure  30.  MDCCA  of  non-intersecting  slices. 


Papers  of  our  novel  methods  have  been  published  in  the  Proceedings  of  the  2006  IEEE  International 
Symposium  for  Biomedical  Imaging:  From  Nana  to  Macro  and  the  Proceedings  of  the  IEEE 
International  Conference  on  Image  Processing  2006.  A  journal  papers  of  our  despeckling,  surface 
segmentation,  and  volume  calculation  (MDCCA)  methods  have  been  submitted  to  Elsevier  DSP, 

Elsevier  Computerized  Medical  Imaging  and  Graphics,  and  IEEE  Transactions  in  Medical  Imaging 

Excellent  progress  with  respect  to  Aim  3  has  been  made  in  the  third  year  of  the  grant. 

The  fourth  Aim  from  the  Statement  of  Work  relates  to  a  small  scale  clinical  validation  in  collaboration 
with  partially  funded  University  of  Virginia  collaborator  -  Dr.  Dan  Theodorescu.  In  the  third  year  of  the 
grant  we  worked  to  resolve  issues  with  obtaining  a  joint  US  DoD  and  UVA  IRB  approval.  We  believe  that 
the  protocol  and  patient  consent  forms  are  in  order  and  are  pending  acceptance.  We  filed  for  a  no  cost 
extension  to  allow  for  completion  of  the  proposed  research  -  i.e.  the  small  scale  clinical  validation.  The 
small  scale  clinical  validation  and  any  subsequent  report/  paper  will  probably  mark  the  conclusion  of  the 
funded  and  reportable  work  in  this  grant. 


33 


KEY  RESEARCH  ACCOMPLISHMENTS 

•  We  have  largely  completed  Aims  1 ,  2  and  3.  Specifically,  we  have  designed  and  had  fabricated  a 
very  high  resolution  transrectal  ultrasound  transducer  array  for  high  resolution  prostate  imaging. 

•  We  have  integrated  the  new  transducer  with  an  ultrasound  scanner  and  an  automated  injection 
stage  to  realize  an  accurate  elastographic  imaging  device. 

•  We  have  tested  the  3D  and  elastographic  imaging  capability  of  the  transducer  /  scanner. 

•  We  have  developed  new  ultrasound  speckle  reduction  techniques  suitable  for  prostate  ultrasound 

•  We  have  developed  new  ultrasound  signal  tracking  techniques  that  improve  the  quality  of 
ultrasound  derived  3D  measures  of  elastic  anomaly  as  might  be  very  useful  in  the  more  sensitive  and 
more  specific  detection  of  prostate  cancer. 
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CONCLUSIONS 

Our  prostate  imaging  approach  combines  using  an  I-Beam  transducer  with  3D  capability, 
elasticity  imaging  and  test  on  a  prototype  using  a  prostate  tissue-mimicking  phantom.  The 
prostate  strain  imaging  performed  here  using  a  slightly  inflated  sheath  over  the  transrectal 
transducer  significantly  enhanced  tumor  visibility  (a  hard  inclusion  in  the  phantom).  (The  lesion 
was  nearly  invisible  in  the  regular  B-mode  image.)  The  I-Beam  transducer  enabled 
reconstruction  of  discrete  2D  image  acquisitions  into  regular  3D  grid  space,  and  thus  the 
tumor  was  rendered  in  3D.  The  volume  calculated  for  this  tumor  had  an  error  of  approximately 
11%  compared  to  the  actual  (independently  determined)  volume.  Additionally,  we  have  made 
significant  progress  in  the  area  of  image  pre-processing  (i.e.  speckle  reduction)  and  in  image 
quantification  (i.e.  feature  segmentation).  These  image  processing  contributions  significantly 
enhance  the  practical  utility  of  our  technique  since  they  hold  the  promise  of  accelerating  the 
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prostate  cancer  diagnostic  task  and  reducing  intra  and  inter  operator  variability.  Reducing 
variability  has  significance  since  serial  analysis  of  cancer  growth  or  remission  is  dependent  on 
accurate  and  repeatable  measures  of  prostate  volume.  Since  we  are  able  to  measure  volumes 
directly,  rather  than  extrapolating  volume  from  a  length  dimension  or  cross-sectional  area,  our 
image  contributions  are  well-matched  and  complement  our  contributions  in  3D  and 
elastographic  imaging. 

“So  what?  Section” 

Our  improved  prostate  ultrasound  imaging  techniques  have  the  potential  to  detect  prostate 
cancer  earlier  and  with  more  reliability  (i.e.  improved  sensitivity  and  specificity).  In  this  way,  we 
believe  that  our  research  has  a  useful  public  health  contribution. 
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